Transitions between Inherent Structures in Water 
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The energy landscape approach has been useful to help 
understand the dynamic properties of supercooled liquids and 
the connection between these properties and thermodynam- 
ics. The analysis in numerical models of the inherent struc- 
ture (IS) trajectories — the set of local minima visited by 
the liquid — offers the possibility of filtering out the vibra- 
tional component of the motion of the system on the poten- 
tial energy surface and thereby resolving the slow structural 
component more efficiently. Here we report an analysis of an 
IS trajectory for a widely-studied water model, focusing on 
the changes in hydrogen bond connectivity that give rise to 
many IS separated by relatively small energy barriers. We 
find that while the system travels through these IS, the struc- 
ture of the bond network continuously modifies, exchanging 
linear bonds for bifurcated bonds and usually reversing the ex- 
change to return to nearly the same initial configuration. For 
the 216 molecule system we investigate, the time scale of these 
transitions is as small as the simulation time scale (« 1 fs). 
Hence for water, the transitions between each of these IS is 
relatively small and eventual relaxation of the system occurs 
only by many of these transitions. We find that during IS 
changes, the molecules with the greatest displacements move 
in small "clusters" of 1-10 molecules with displacements of 
~ 0.02 — 0.2 nm, not unlike simpler liquids. However, for 
water these clusters appear to be somewhat more branched 
than the linear "string-like" clusters formed in a supercooled 
Lennard- Jones system found by Glotzer and her collabora- 
tors. 



I. INTRODUCTION 

It is believed that the properties of liquids can be un- 
derstood as motion of the system in a high-dimensional 
complex potential energy surface (PES) As a liquid 

is cooled toward the glassy state, the system is increas- 
ingly found near local potential energy minima, called 
inherent structure (IS) configurations 0. As tempera- 
ture decreases, the description of the dynamics in terms 
of motion on the PES becomes increasingly appropriate. 
In this description, in the glassy state, the system is lo- 
calized in one of the potential energy basins [fzHTT 



While such a picture of liquid dynamics is difficult to 
verify experimentally, computer simulation offers an ex- 
cellent opportunity to explore these ideas. For a pre- 
defined liquid potential, a liquid trajectory can be gen- 
erated via molecular dynamics simulation and the lo- 
cal potential energy minima can be evaluated by an en- 
ergy minimization method B. With this procedure, the 
motion in phase space is converted into a minimum-to- 
minimum trajectory, or IS trajectory. A general picture 
of the system moving among a set of basins surrounding 
the multitude of local minima has evolved. More specifi- 
cally, simulations have shown that both the depth of the 
minima sampled by the system, as well as the number 
of these minima, decrease on cooling [|[ |l^, [ll]] . Simula- 
tions have also shown that below a crossover temperature 
T x , only rare fluctuations bring the system to a saddle 
point and hence activated processes become important 
for relaxation of the liquid fipj, |l2|- |Ii|[ . For the system 
we study, T x coincides numerically [|| |ld|-|l7j with the 
critical temperature identified by mode coupling theory 
(MCT), which has been widely used to understand the 
dynamics of liquids on "weak" supercooling (the T range 
where characteristic relaxation times approach ps 10~ 8 s) 
& 

The description of the real motion of the system as 
an IS trajectory becomes a powerful way of separating 
the vibrational contribution, responsible for the thermal 
broadening of instantaneous measurements from the slow 
structural component |l9| . Such an approach becomes 
even more powerful below T x , since most of the instan- 
taneous configurations are far from saddles, making cor- 
relation functions calculated from the IS trajectory fully 
account for the a-relaxation dynamics M. 

Here we study the IS trajectory below T x for the ex- 
tended simple point charge (SPC/E) potential^)]], a well- 
studied model for water. The dynamics of the SPC/E 
model have been shown to be consistent with the predic- 
tions of MCT Additionally, the PES of SPC/E 
has been studied in detail above T x and a thermody- 
namic description of the supercooled states based on the 
PES has been presented|l^,|3). Here we focus on the ge- 
ometrical properties of the motion, once the vibrational 
component is subtracted. The possibility of performing 
such a study below T x , with a very fine time coarse grain- 
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ing, allows us to examine the structural changes that 
accompany the basin transitions and to describe an ele- 
mentary step of the diffusive process in terms of hydrogen 
bond network rearrangement. 

This work is organized as follows. In Sec. |j] we pro- 
vide the simulation details. In Sec. Ill we analyze the 



characterized by an energy change of w 10 — 20 kJ/mol 
and an oxygen atom square displacement of the order of 
0.01 A 2 ; they appear to constitute the elementary step 
underlying the diffusional process in the system. 

We note that it is impossible to identify the transi- 
tion unless the quenching time step is of the order of 
the simulation time step. This is in distinct contrast 
to a Lennard-Jones liquid, where such small continuous 
changes are not present [Q. The difference, we will show, 
is attributable to the hydrogen bonds. This time scale for 
IS transitions is in accord with other work on the TIPS2 
water model, at T = 298 K @. 

In Fig. ||, we see the sharp changes in Eis(t) coincide 
with the sharp changes in (r 2 (t)). This confirms that the 
Our results are based on molecular dynamics Simula- system is repeatedly visiting specific configurations, since 



displacement of the molecules in these IS-transitions and 
in Sec. IV we study the corresponding changes in the 
hydrogen-bond (HB) network. Finally, in Sec. ^ we 
present a brief summary. 



II. SIMULATION 



tions of the SPC/E model pj] of water for 216 molecules, 
at fixed density p — 1 g/cnr . The numerical procedure 
is described in Ref.^2|. The integration time step St 
is 2 fs. The mode coupling temperature for this den- 
sity is Tmct ~ 194K" p^] . We analyze trajectories at 
T = 180 K, so that the system is in the deep supercooled 
liquid state. At this temperature, the diffusion coeffi- 
cient is more than four orders of magnitude smaller than 
its value at T = 300 K and only a few molecules move 
significantly (with displacements larger than 0.025 nm) 
at each simulation time stepj24j. Our system is started 
from equilibrated configurations at 190 K, which relax 
for nearly 920 ns at 180 K before we record and analyze 
the trajectory. At such a low temperature, a slow aging 
in the trajectory could be present, however, the aging 
should not affect the qualitative picture we present. 

We have generated one trajectory of 30 ns, sampling 
configurations at each 1 ps. For each configuration we 
find the corresponding IS using conjugate gradient min- 
imization. In this way, we obtain 30, 000 configurations 
with the corresponding IS. Since we could miss some IS 
transitions with 1 ps sampling, we also ran four indepen- 
dent 20 ps simulations sampling the IS at 4 fs. In this 
way, we obtain another 20,000 configurations with the 
corresponding IS. 



III. INHERENT STRUCTURE TRAJECTORIES 

Figure ^ shows an example of the potential energy 
of the IS, Eis(t), and the mean square displacement of 
the oxygen atoms (r 2 (t)) starting from a single arbitrary 
starting time. At T = 180 K, the slowest collective re- 
laxation time r a > 200 ns p2| . The IS trajectory in 
Fig. |l| has a mesh of 1 ps and covers a total time of 
30 ns. In this time interval, (r 2 (t)) is about lA 2 , i.e. 
much less than the corresponding value of the average 
nearest neighbor distance of 2.8A. Figure shows an en- 
largement of the IS trajectory using a much smaller time 
mesh (4 fs, two times the simulation time step). Fig. ^ 
shows that changes occur via discrete transitions, with 
an average duration of f» 0.2 ps. The transitions are 



(r (t)) and Eis(t) take on discrete values. The results 
shown in Fig. g imply that the system often returns to 
the original basin because of both the difference in energy 
and the displacement approach zero at the end of the time 
interval. 

To aid in understanding the distribution of the dis- 
placements during the IS changes (such as those between 
the two IS labeled A and B in Fig. |(b)), Fig. | shows 
the displacements u of all 216 individual molecules from 
IS- A to IS-B. We see that there is a relatively small set of 
molecules with a large displacement. A snapshot of the 
eight molecules with the largest displacement is shown in 
Fig. ^. Interestingly, we find that this set of molecules 
forms a cluster of bonded molecules. Indeed, for all cases 
studied, we found that the set of molecules which dis- 
place most form a cluster of bonded molecules. The ob- 
served clustering phenomenon characterizes the IS tran- 
sitions in water and can be interpreted as the analog of 
the string-like motion observed in simple atomistic liq- 
uids and connected to the presence of dynamical het- 
erogeneities j25|. Similar results were found by Ohmine 
et al. using the TIP4P and TIPS2 models for water 

To characterize the distribution of individual molecu- 
lar displacements between different IS more carefully, we 
show the distribution of displacements u of the oxygen 
atoms P(u, 5t — 4 fs) in Fig. 0. The time difference St 
is two times the integration time step^Tj. We collected 
data from all transitions in four independent 10 ps trajec- 
tories for the results reported in Fig. |j| Note that P{u) 
was previously studied by Schr0der et al. for a binary 
Lennard-Jones (LJ) mixture [Q. 

We find that the distribution of displacements is mono- 
tonic and does not allow us to detect a characteristic 
length which could help in distinguishing between diffu- 
sive and non-diffusive basin changes. P(u,t) for t = 5t 
shows an apparent power-law with a negative slope of 
about —2.8 ± 0.2, followed by an exponential tail. To 
highlight the exponential tail, we present the data in a 
linear- log plot Fig. || (upper panel), showing that the tail 
in P(u, 5t — 4 fs) is mostly due to these highly "mobile 
molecules". The characteristic length of the exponen- 
tial tail is 0.2A, about 15 times smaller than the nearest 
neighbor oxygen distance of 2.8A and in agreement with 
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the finding of Schr0der et al. for the LJ case, where the 
characteristic length was about 1 /8 of the nearest neigh- 
bor distance. 

Schr0der et al. interpreted the power-law contribu- 
tion to P(u, St) as a response of the system to the diffus- 
ing molecules. In elasticity theory [f28|, a local displace- 
ment at the origin produces a distribution of displace- 
ments which scales with the distance R from the origin 
as R~ 2 . Hence the distribution of displacements u scales 
as u~ 2 ' 5 ^8|. In contrast with the LJ case, the exponent 
we find is slightly larger than 2.5, a discrepancy which 
may be related to the application of a continuum theory 
at a length scale where molecular details are still rele- 
vant. In the spirit of elasticity theory, the deviation of 
the distribution at small u values from the power law 
arises from the cutoff introduced by the finite size of the 
simulation. Indeed, small u values are produced at large 
R and hence are missing in a finite simulation box. 

Although a subset of "highly mobile" molecules is iden- 
tifiable using a pre-defined threshold value in a single 
basin change, there is no unambiguous general criterion 
for identifying the molecules responsible for a single basin 
transition. 

An important open question is the relation between 
the displacement distribution functions for t — St (Fig. ||) 
and the same distribution evaluated for an arbitrary time 
interval t. A full knowledge of such relation may shed 
light on the elementary stochastic process which better 
describes the dynamics of the deep supercooled state. 
This calculation requires a significant computational in- 
vestment due to the very large number of IS configura- 
tions which need to be evaluated. Despite such difficul- 
ties, it is a promising research line for the near future. 

IV. HB NETWORK RESTRUCTURING 

Liquid water is an interesting system for studying the 
physical processes that accompany basin transitions. The 
unusual dynamic and thermodynamic properties of water 
are believed to be connected to the microscopic behav- 
ior of hydrogen bonding. Many experiments suggest that 
this tetrahedral HB network has defects, such as an ex- 
tra (fifth) molecule in the first coordination shell ]29|, |30|. 
Indeed, such 5-coordinated molecules have been directly 
identified in simulations and the defects were found to be 
a catalyst for motion in the system[^l[, making an ob- 
vious possible connection between network defects pro- 
moting diffusion and the basin transitions that give rise 
to diffusive motion of that system. 

To obtain a physical picture of the IS transitions for 
water and to hopefully better understand the source of 
these transitions, we focus on the small changes in the 
HB network along a minimum to minimum trajectory. 
Analysis of the HB network based on IS configurations 
has shown that local PES minima are characterized by 
both linear bonds (LB) and bifurcated bonds (BB) pi] ] 



whose fraction is both temperature and density depen- 
dent. The HB network tends toward a perfect random 
tetrahedral network on cooling, or on lowering the pres- 
sure [E3[. Here we emphasize the changes in the network 
associated with IS-transitions. 

To quantify the bonding changes, it is necessary to 
employ an arbitrary bond definition. Previous works 
have indicated that several definitions provide physically 
reasonable results. Here we use the definition that two 
molecules are bonded if their oxygen-oxygen distance is 
less than 3.5A and their mutual potential energy is nega- 
tive |3l| ]. Using this definition, we obtain the equilibrium 
distribution of the potential energy for the LB and BB, 
shown in Fig. |^. 

These results agree with previous findings, based on 
different models p3-p4]. The distribution of BB ener- 
gies is bimodal with peaks at roughly —6 kJ/mol or 
—22.5 kJ/mol, while the LB energy distribution is uni- 
modal with a peak at roughly —24 kJ/mol. Therefore, 
the energy associated with a change in the HB network 
due to losing one LB and creating two BB ranges roughly 
from —21 kJ/mol to 12 kJ/mol, depending on which of 
the two possible BB are created. The relative intensities 
of the peaks of the BB energy distribution suggest that 
such a mechanism would more likely lead to an increase 
in the overall energy. The comparison of the HB-encrgy 
changes with those found in Fig. [I] suggest that inter- 
changes between LB and BB can explain the existence of 
a multitude of IS separated by very small energy barriers. 
More specifically, we hypothetize that changes increas- 
ing Ejs found in Fig. [l] are due to processes LB— >BB, 
while the changes decreasing Ejs are due to the pro- 
cesses BB— >LB. To confirm or reject this hypothesis, we 
study the distributions of LB and BB for IS just before 
and just after positive and negative "jumps" in the po- 
tential energy (Fig. with energy larger than 9 kJ/ mol 
p5| . These IS are schematically shown in Fig. 0(a). The 
corresponding distributions for LB and BB are shown in 
Fig. @(b) and@(c). We see that during PE increases, the 
average number of LB decreases while the number of BB 
increases; the opposite situation occurs when the poten- 
tial energy decreases. The distribution for BB has peaks 
at even numbers of BB, which we expect since for each 
LB lost, two BB appear, also implying that the distri- 
bution should be zero for odd numbers of BB |3(| . The 
anti-correlation of the BB and LB changes is a strong 
indicator that a mechanism whereby the system accesses 
higher energy states is via LB^BB transitions. 

To reinforce the hypothesis that the basin change is 
associated with a restructuring of the local connectivity, 
Fig. ^ shows the number of molecules with a coordina- 
tion number equal to 3, 4 or 5 as a function of time for a 
characteristic time interval and contrasts this data with 
the time dependence of (r 2 (t)). A clear anti-correlation 
is observed between the time dependence of the number 
of 3 and 5-coordinated molecules compared to the time 
dependence of the 4-coordinated molecules, supporting 
the proposed interchange mechanism. The fact that in- 
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creases in the fraction of BB coincides with changes in 
(r 2 (t)) supports the expectation that motion is a result 
of network imperfections. 

V. SUMMARY 

We have presented a detailed analysis of the motion 
on the PES of a 216 water molecule system interacting 
via the SPC/E potential in a deeply supercooled state, 
below the MCT temperature. At this temperature, the 
system populates basins of the local minima. At these 
conditions, the analysis of the IS trajectories provides a 
very clean description of the slow alpha relaxation pro- 
cess filtering nearly all vibrational motions. 

We have shown that an inherent structure transition 
is observed about every 0.2 ps. It is the collection of 
these numerous small transitions that gives rise to the 
structural relaxation of the system. These fast transi- 
tions are characterized by a broad distribution of indi- 
vidual molecule displacements, without a clear charac- 
teristic length. Future work must address the issue of 
the prediction of P(u, t) from the knowledge of P(u, St). 

We perform an analysis of the geometry of the indi- 
vidual event and find that the most mobile molecules are 
clustered. The analysis of the changes in HB connectiv- 
ity associated with IS changes reveals that these transi- 
tions are associated with the breaking and reformation 
of HB. This result is in accord with work of Ohmine 
and colleagues @|§ for TIP4P. We have shown that 
the transitions associated with an increase in the energy 
correspond to the breaking of linear bonds and to the 
simultaneous formation of bifurcated bonds. Similarly, 
the transitions associated with a decrease in the energy 
correspond to the breaking of bifurcated bonds and to 
the simultaneous formation of linear bonds. 

This result supports the hypothesis that the linear to 
bifurcated transition can be considered as an elementary 
step in the rearrangement of the HB network. 
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FIG. 1. (a) Mean square displacement and (b) IS energy 
for the inherent structures as a function of time for the stud- 
ied 216-molecule system. The time interval between adjacent 
IS in both figures is 1 ps. While it is possible to track IS 
transitions from the potential energy, it is not the case for 
the mean square displacement. Note the amplitude of the 
peaks of the potential energy is ~ 10 — 20 kj/mol, the same 
order of magnitude of a HB energy. 
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FIG. 2. (a) IS energy and (b) mean square displacement for 
the IS obtained using a sampling interval of 4 fs, two times 
the simulation time step. The correlation between Eis and 
(r 2 (t)) is evident. Also, we see that it is necessary to sample 
IS with a mesh of the order of the simulation time step to 
detect all the IS visited by the system. 
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FIG. 3. Displacement of each molecule in the transition 
from IS- A to IS-B shown in Fig. |(b). 




FIG. 4. Snapshot of the system in the IS labeled A in 
Fig. ^(b). Only the eight molecules with displacement larger 
than 0.025 nm are shown here. Hydrogen-bonded molecules 
are connected by tubes. Note that all 8 molecules are nearby 
and form a cluster, which unlike the LJ case, are bonded and 
less string-like. 
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FIG. 6. Distribution function for the pair interaction en- 




ergy Vij for the LB and BB in the IS. Note the bimodal dis- 
tribution for BB and unimodal distribution for LB. 
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FIG. 5. Distribution of displacements u of the oxygen 
atoms between IS changes, P(u,i), evaluated at t = 4 fs. We 
show the distribution of all the displacements of the molecules 
(squares) and the distribution for the largest displacement 
in an IS-transition (circles). For comparison, we divided the 
largest displacement distribution by the number of molecules. 
We present both (a) log-log and (b) linear-log graphs to show 
the power law behavior for u ~ 0.01 nm and the exponen- 
tial tail of the distributions. We also show the prediction of 
elasticity theory P(u) ~ u~ 2 ' 5 . 
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FIG. 7. Panel (a) shows schematically the separation of IS 
used to meassure the probability distributions for (b) LB and 
(c) BB. Circles correspond to the distribution calculated over 
IS just before an increase in potential energy larger than 9 
kj/mol, i.e. when the system is in an IS A shown in (a). 
Squares correspond to the distribution calculated over IS just 
after an increase in potential energy larger than 9 kj/mol, i.e. 
when the system is in an IS B shown in (a). Similarly, dia- 
monds correspond to the distribution calculated over IS just 
before a decrease in potential energy larger than 9 kj/mol, 
when the system is in an IS C shown in (a). Triangles cor- 
respond to the distribution calculated over IS just after a 
decrease in potential energy larger than 9 kj/mol, when the 
system is in an IS D shown in (a). 
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FIG. 8. Number of molecules with coordination number 
CN equals 3,4 and 5 versus time, for the IS corresponding to 
Fig. |. The plot for CN = 4 is shifted down by 400 units 
for better comparison. Also shown is (r 2 (t)}. We see how the 
tetrahedral network acquires both types of defects (CN equals 
to 3 and 5) while the system explores different IS. 



